3 Algoritmo

set.seed(288) ## setando a semente para gerar sempre os mesmos resultados da amostra
amostra <- rgamma(112,420,10) # gerando amostra da gamma

A amostra segue uma distribuição \(gama(\alpha = 420,\beta = 10)\), logo esses são os verdadeiros parâmetros que tentaremos estimar através do método bayesiano, note que \(E(Y) = \frac{\alpha}{\beta} = 42\), logo nesse exemplo criado a empresa em média não cumpre com o plano contratado de 50 Megabits. Obs: Para a análise vamos considerar todos os parâmetros desconhecidos e seus valores reais não serão mencionados a partir daqui.

n = length(amostra) # tamanho da amostra
media = mean(amostra) # media amostral
variancia = var(amostra) # variância amostral
soma_y = sum(amostra) # soma da amostra
soma_log_y <- sum(log(amostra)) # soma do log da amostra
media ## valor da media
## [1] 42.09928

Aqui foi calculado algumas estatísticas que serão utéis para mais adiante, é possivel notar que a média foi bem consistente com seu verdadeiro parâmetro, indicando uma boa qualidade da amostra.

3.1 Escolha das Prioris

Como o objetivo é averiguar se a operadora está ofertando uma conexão conforme o previsto, se escolheu prioris para alpha e beta de tal maneira que a operadora não possa atribuir o resultado negatativo do teste às prioris. Abaixo foi mostrado esse Gráfico de uma \(gama(\alpha = 50,\beta = 1)\) e perguntado se representa bem a variabilidade da internet, João concordou com a distribuição. A distribuição segue a premissa de que a operadora distribui em média uma velocidade de conexão mínima, portanto com esperança igual a 50, logo a operadora não pode reclamar de viés da priori.

alpha <- 50
beta <- 1
x <- seq(qgamma(0.0001,shape = alpha,rate = beta),qgamma(0.9999,shape = alpha,rate = beta),0.01)
p <- plot_ly(x = x, y = dgamma(x,shape = alpha, rate = beta), type = 'scatter', mode = 'lines', fill = 'tozeroy') 
p
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

Logos as prioris \(\alpha_0\) e \(\beta_0\) foram escolhidas de tal maneira que \(E(\alpha_0) = 10\) e \(E(\beta_0) = 1/5\) e com variâncias bem altas.

# constante para alpha0, note que quando menor k1 maior será a variancia!
k1=1/5

## Priori para alpha
lambda_alpha = 50*k1 # lambda da priori alpha
nu_alpha = k1 # nu da priori alpha


x <- seq(0,qgamma(0.995,shape = lambda_alpha,rate = nu_alpha),0.01)
p <- plot_ly(x = x, y = dgamma(x,shape = lambda_alpha, rate = nu_alpha), type = 'scatter', mode = 'lines', fill = 'tozeroy') 
p
## Priori para beta
k2 = 1/2
lambda_beta = k2 # lambda da priori beta
nu_beta = k2 # nu da priori beta
x <- seq(0,qgamma(0.99,shape = lambda_beta,rate = nu_beta),0.0001)
p <- plot_ly(x = x, y = dgamma(x,shape = lambda_beta, rate = nu_beta), type = 'scatter', mode = 'lines', fill = 'tozeroy') 
p

Pelos gráficos conseguimos perceber que as prioris são pouco informativas como queriamos.

### Passo 1
n_simu <- 300000 # número de simulações
simu <- matrix(nrow = n_simu,ncol = 3) # inicializando matriz com as simulações
colnames(simu) <- c("alpha","beta","x") # nome das colunas da matriz
#simu[1,1:3] <- c(media^2/variancia,media/variancia,media) # calculando a primeira linha
simu[1,1:3] <- c(1000,100,media) # calculando a primeira linha
var2 = 10 ## variancia da gama que servira como distribuicao proposta

simular_proposta <- function(mu){  ## gera valor de uma gamma com media igual a mu e variancia igual a var2
  rgamma(1,shape = mu^2/var2,rate = mu/var2)
}

proposta <- function(y,mu){ ## calcula a densidade da distribuicao proposta no ponto y dado mu
  dgamma(y,shape = mu^2/var2,rate = mu/var2)
}
priori_alpha <- function(y){ ## calcula a densidade para a priori de alpha
  dgamma(y,lambda_alpha,nu_alpha)  + (dgamma(y,lambda_alpha,nu_alpha) == 0)*.000001 ## esse ultimo termo é só para valores de priori igual a 0 não dar problema com log
}

log_vero <- function(alpha,beta){ ## log verossimilhança 
  n*alpha*log(beta) - n*lgamma(alpha) + (alpha-1)*soma_log_y
}
## só para calcular o indice de aceitação de alpha
aceita = c()

set.seed(2465)

for (i in 2:n_simu) {
  ## pegando os valores anteriores da cadeia
  alpha_anterior <- simu[i-1,1] 
  beta_anterior <- simu[i-1,2]
  
  ### Metropolis Hasting para gerar alpha
  
  alpha_proposto <- simular_proposta(mu = alpha_anterior) ## simulando um valor proposto para alpha da distribuicao gama com media igual ao alpha do passo anterior
  
  ## calculo do da log_verossimilhanca no ponto (alpha proposto,beta do passo anterior)
  loglik.atual <- log_vero(alpha_proposto,beta_anterior) 
  ## calculo do da log_verossimilhanca no ponto (alpha passo anterior,beta do passo anterior)
  loglik.anterior <- log_vero(alpha_anterior,beta_anterior)
  
  ## calculo da densidade da log priori no ponto alpha proposto/alpha anterior
  log.prior.alpha.atual <- log(priori_alpha(alpha_proposto))
  log.prior.alpha.anterior <- log(priori_alpha(alpha_anterior))
  
  ## calculo da densidade da log da distribuição proposta no ponto do alpha proposto dado o alpha anterior e do alpha anterior dado o alpha proposto
  log.proposta.atual <- log(proposta(alpha_proposto,alpha_anterior))
  log.proposta.anterior <- log(proposta(alpha_anterior,alpha_proposto))
  
  logpi.atual <- loglik.atual + log.prior.alpha.atual - log.proposta.atual
  logpi.anterior <- loglik.anterior + log.prior.alpha.anterior - log.proposta.anterior
   
  u <- runif(1)
  aceita[i-1] <- log(u) <= (logpi.atual - logpi.anterior)
  
  alpha_atual <- aceita[i-1]*alpha_proposto + (1- aceita[i-1])*alpha_anterior
  
  ###############
  # Atualizar Beta
  ###############
  
  ## gerando o valor de beta pela posteriori condicional
  beta_atual <- rgamma(1, alpha_atual*n + lambda_beta, rate = soma_y + nu_beta)  
  
  ## gerando um valor de gamma (desnecessário, mas quis deixar :P)
  x_atual <- rgamma(1, alpha_atual,beta_atual)
  
  ## Guardando os valores 
  simu[i,] <- c(alpha_atual, beta_atual,x_atual)
  
}

Taxa de aceitação da cadeia

taxa_aceitacao <- sum(aceita)/length(aceita)
taxa_aceitacao
## [1] 0.4533715

3.2 Resultado da cadeia

cadeia1 <- simu[,1:2]

theta = as.mcmc(cadeia1)
plot(theta)

Com esse primeiro gráfico podemos notar que a cadeia convergiu, apesar que tenha demorado um pouco.

theta2=window(as.mcmc(theta),start=5000,thin=1)
plot(theta2)

Retirando as 5000 primeiras observações da cadeia podemos perceber que não é mais necessário retirar mais observações

autocorr.plot(theta2,lag.max = 5000)

No gráfico de autocorrelação podemos perceber um grande, as observações são muito correlacionados e só começam a diminuir com um número muito grande de espaçamento entre as observações, com aproximadamente 1500 observações de espaçamento as observações se tornam independentes, o problema é que com isso precisamos de 1500 simulações para gerar um valor da posteriori independente, que causa um enormo custo computacional.

ps: Nos tentamos de muitas maneiras de diminuir esse problema, mas nenhuma foi efetiva :(

theta2=window(as.mcmc(theta),start=5000,thin=1500)
plot(theta2)

autocorr.plot(theta2,lag.max = 10)

3.3 Interpretação dos Resultados

summary(theta2)
## 
## Iterations = 5000:299000
## Thinning interval = 1500 
## Number of chains = 1 
## Sample size per chain = 197 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##          Mean      SD Naive SE Time-series SE
## alpha 210.735 25.2940  1.80212        1.80212
## beta    5.004  0.6041  0.04304        0.04304
## 
## 2. Quantiles for each variable:
## 
##          2.5%    25%     50%     75%   97.5%
## alpha 163.569 192.89 210.279 228.851 258.058
## beta    3.868   4.59   4.983   5.429   6.123

A esperança para \(\alpha\) tem valor estimado igual a 210.735 e mediana igual a 210.279, o intervalo de credibilidade central de 0.95 de probabilidade é dado por (163.569;258.058). Para \(\beta\) a esperança é igual a 5.004 e mediana igual a 4.983 com intervalo de credibilidade central de 0.95 de probabilidade dado por (3.868;6.123).

Se considerarmos a esperanças dos parâmetros então \(E(Y)= \frac{\alpha}{\beta} = 42.11331\), como esse valor envolve dois parametros é dificil gerar um valor para essa estimativa, talvez fosse o caso de simular com uma noma parametrização.

3.4 Gerar com outra distribuição proposta para \(\alpha\)

### Passo 1
n_simu <- 300000 # número de simulações
simu2 <- matrix(nrow = n_simu,ncol = 3) # inicializando matriz com as simulações
colnames(simu2) <- c("alpha","beta","x") # nome das colunas da matriz
#simu[1,1:3] <- c(media^2/variancia,media/variancia,media) # calculando a primeira linha
simu2[1,1:3] <- c(1000,100,media) # calculando a primeira linha
simular_proposta2 <- function(mu){  ## gera valor de uma qui quadrado com media igual a mu
  rchisq(1,df=mu)
}

proposta2 <- function(y,mu){ ## calcula a densidade da distribuicao proposta no ponto y dado mu
  dchisq(y,df = mu)
}
## só para calcular o indice de aceitação de alpha
aceita = c()

set.seed(2465)

for (i in 2:n_simu) {
  ## pegando os valores anteriores da cadeia
  alpha_anterior <- simu2[i-1,1] 
  beta_anterior <- simu2[i-1,2]
  
  ### Metropolis Hasting para gerar alpha
  
  alpha_proposto <- simular_proposta(mu = alpha_anterior) ## simu2lando um valor proposto para alpha da distribuicao gama com media igual ao alpha do passo anterior
  
  ## calculo do da log_verossimilhanca no ponto (alpha proposto,beta do passo anterior)
  loglik.atual <- log_vero(alpha_proposto,beta_anterior) 
  ## calculo do da log_verossimilhanca no ponto (alpha passo anterior,beta do passo anterior)
  loglik.anterior <- log_vero(alpha_anterior,beta_anterior)
  
  ## calculo da densidade da log priori no ponto alpha proposto/alpha anterior
  log.prior.alpha.atual <- log(priori_alpha(alpha_proposto))
  log.prior.alpha.anterior <- log(priori_alpha(alpha_anterior))
  
  ## calculo da densidade da log da distribuição proposta no ponto do alpha proposto dado o alpha anterior e do alpha anterior dado o alpha proposto
  log.proposta.atual <- log(proposta(alpha_proposto,alpha_anterior))
  log.proposta.anterior <- log(proposta(alpha_anterior,alpha_proposto))
  
  logpi.atual <- loglik.atual + log.prior.alpha.atual - log.proposta.atual
  logpi.anterior <- loglik.anterior + log.prior.alpha.anterior - log.proposta.anterior
   
  u <- runif(1)
  aceita[i-1] <- log(u) <= (logpi.atual - logpi.anterior)
  
  alpha_atual <- aceita[i-1]*alpha_proposto + (1- aceita[i-1])*alpha_anterior
  
  ###############
  # Atualizar Beta
  ###############
  
  ## gerando o valor de beta pela posteriori condicional
  beta_atual <- rgamma(1, alpha_atual*n + lambda_beta, rate = soma_y + nu_beta)  
  
  ## gerando um valor de gamma (desnecessário, mas quis deixar :P)
  x_atual <- rgamma(1, alpha_atual,beta_atual)
  
  ## Guardando os valores 
  simu2[i,] <- c(alpha_atual, beta_atual,x_atual)
  
}
cadeia2 <- simu2[,1:2]

theta3 = as.mcmc(cadeia2)
plot(theta3)

autocorr.plot(theta3,lag.max = 5000)

theta4=window(as.mcmc(theta3),start=5000,thin=3500)
plot(theta4)

autocorr.plot(theta4,lag.max = 10)

summary(theta4)
## 
## Iterations = 5000:299000
## Thinning interval = 3500 
## Number of chains = 1 
## Sample size per chain = 85 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##         Mean      SD Naive SE Time-series SE
## alpha 213.54 24.4106  2.64770        2.64770
## beta    5.07  0.5869  0.06365        0.06365
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%     50%     75%   97.5%
## alpha 172.006 194.555 210.408 228.851 261.990
## beta    4.085   4.607   4.999   5.446   6.229

Nesse o gráfico de autocorrelacao ficou pior ainda, precisando de um espaçamento de 3500 observações.